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Abstract 

We give a thorough Biorthonormal Matrix-Product-State (BMPS) analysis of the 
Transfer-Matrix Renormalization-Group (TMRG) for non-Hermitian matrices in the 
thermodynamic limit. The BMPS is built on a dual series of reduced biorthonormal 
bases for the left and right Perron states of a non-Hermitian matrix. We propose 
two alternative infinite-size Biorthonormal TMRG (iBTMRG) algorithms and compare 
their numerical performance in both finite and infinite systems. We show that both 
iBTMRGs produce a dual infinite-BMPS (iBMPS) which are translationally invariant in 
the thermodynamic limit. We also develop an efficient wave function transformation of 
the iBTMRG, an analogy of McCulloch in the infinite-DMRG arXiv:0804.2509 (2008)], 
to predict the wave function as the lattice size is increased. The resulting iBMPS 
allows for probing bulk properties of the system in the thermodynamic limit without 
boundary effects and allows for reducing the computational cost to be independent of 
the lattice size, which are illustrated by calculating the magnetization as a function of 
the temperature and the critical spin-spin correlation in the thermodynamic limit for a 
2D classical Ising model. 

Keywords: Transfer-matrix renormalization-group. Density-matrix renormalization- 
group, Biorthonormal matrix-product-state. Correlation function. 
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1 Introduction 

It has been widely understood that the Matrix-Product-State (MPS) ansatz constitutes the 
basis of many numerical algorithms in computational physics, notably the Density-Matrix 
Renormalization-Group (DMRG) [1] and the Time-Evolving Block Decimation (TEBD) [2]. 
There are two main variants of the DMRG algorithm. The finite-size algorithm was by now 
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realized to compute variationally the ground state of a strongly correlated ID quantum 
system within the class of MPS [3]. While the infinite-size algorithm grows the system 
by adding iteratively one or more sites and reaches asymptotically, at the fixed point, a 
translationally invariant MPS state (invariant under translations of some fixed number of 
lattice sites), it was originally conceived as only for obtaining the initial MPS wave function 
for the finite-size algorithm. Until recently, interest in directly obtaining the translationally 
invariant MPS wave function in the thermodynamic limit was rekindled by the ideas of 
infinite-TEBD (iTEBD) of Vidal [4] and infinite-DMRG (iDMRG) of McCulloch [5]. Both 
approaches allow for probing bulk properties of the system in the thermodynamic limit 
without the influence of boundary conditions and the resulting infinite-MPS (IMPS) can 
reduce the computational cost of real-time simulations to be independent of the lattice size. 
The main difference between these two methods lies in the scheme for the local update of 
the tensors in the MPS. The iTEBD applies a single bond evolution operator to each site 
simultaneously which amounts to a power method while the iDMRG obtains the center 
matrix variationally by using a very efficient local eigensolver. 

On the other hand, after nearly two decades of development, the application of DMRG 
methods has spread over a great variety of fields. One major branch of the DMRG is 
the Transfer-Matrix Renormalization-Group (TMRG) which can be applied to, e.g., the 
strongly correlated classical systems [6], the thermodynamics of ID quantum systems at 
finite temperature [7], the stochastic transfer matrix of a cellular automaton [8], the non- 
equilibrium systems in statistical physics [9] , and the general Markov random field in image 
modeling [10]. Different from the DMRG, the TMRG usually deals with non-Hermitian 
matrices which involve much more numerical demand in view of the existence of distinct 
left and right Perron states (i.e., the eigenstates associated with the maximum eigenvalue 
which we will refer to as the Perron root) of the transfer matrix. This leads to the selection 
of the reduced density matrix becoming ambiguous. Enss and SchoUwock [8] had provided 
a comparative discussion on several choices of the density matrix proposed in the litera- 
ture for the non-Hermitian TMRG. Unfortunately, all these selections of density matrices 
cause the conventional TMRG fail to fit the framework of MPS analysis. In a recent paper 
[11], the author proposed a new TMRG algorithm called Biorthonormal Transfer-Matrix 
Renormalization-Group (BTMRG) which employed a dual set of biorthonormal bases to 
construct the renormalized transfer matrix and reduced the numerical complexity for non- 
Hermitian matrices to be the same as the Hermitian case. Numerical simulations for a real 
non-Hermitian matrix showed that the BTMRG exhibits significant improvement on the ef- 
ficiency and accuracy than conventional TMRG. Here, a dual biorthonormal bases indicate 
any two sets of vectors {|a)}a=i,--- ,m and ... ^ that satisfy {a\ /3) = 5ap- In the pre- 

vious paper, the BTMRG was described in a traditional formulation. In this paper, we will 
show that the employment of the biorthonormal bases enables the BTMRG to be perfectly 
reformulated within the framework of MPS analysis where the two Perron states can be 
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represented as a dual Biorthonormal Matrix-Product-States (BMPS). We will propose two 
alternative methods for building the BMPS. At the same time of doing so, two natural ques- 
tions arise: Does the infinite-size variant of the BTMRG (iBTMRG) create asymptotically 
a dual infinite biorthonormal MPS (iBMPS) that are translationally invariant? Does there 
exist an efficient transformation in iBTMRG, just as the one in iDMRG [5], for predicting 
the wave function as the lattice size is increased? In this paper, both questions will be 
addressed confirmatively. Here, instead of the two-site iDMRG where the non-zero trunca- 
tion of the wave function often leads to undesirable effects of the convergence of the IMPS 
[5], we will use the "single-site" iBTMRG scheme which achieves zero truncation similar 
to the single-site iDMRG. In addition, by using a special E»S»E configuration, the critical 
two-point correlation functions in the thermodynamic limit, especially when the distance 
between the two points approaches infinity, will be shown to be easily obtained. 

The rest of this paper is organized as follows. In section 2, we present a brief review of 
the MPS analysis of the standard DMRG in finite systems and reformulate the BTMRG 
algorithm with the E«S«E configuration proposed in [11] entirely in terms of the BMPS 
language. In section 3, we propose two alternative methods to build the BMPS of the left 
and right Perron state for non-Hermitian matrices and develop an iBTMRG algorithm to 
obtain the asymptotics of the dual iBMPS wave functions that arc translationally invariant. 
We also develop an efficient wave function transformation in iBTMRG for predicting the 
BMPS as the lattice size is increased. Section 4 compares the numerical performance of 
our iBTMRG algorithms in finite systems for a real non-Hermitian matrix of an anisotropic 
Ising model. The magnetization as a function of the temperature and the critical two-point 
correlation function of an isotropic Ising model in the thermodynamic limit are also plotted. 
Finally, in section 5, some conclusions are drawn. 

2 Biorthonormal MPS analysis of BTMRG in finite systems 

The connection between DMRG and MPS was first found by Ostlund and Rommer [12] who 
identified the thermodynamic limit of DMRG with a position-independent matrix product 
wave function. The discovery of the MPS on which the DMRG operates has placed the 
algorithm on a firm footing, provided a deeper understanding, and allowed a concrete and 
easy to manipulate description of the DMRG. The standard DMRG and MPS formulation 
was established on a series of reduced orthonormal bases. However, when considering both 
left and right Perron states of non-Hermitian transfer matrices, the MPS must be built 
on a dual series of reduced biorthonormal bases which we will refer to as biorthonormal 
matrix-product states (BMPS). In this section, we briefiy review the finite-size BTMRG 
algorithm proposed in [11] but reformulate it entirely in terms of a BMPS language. For 
further information about MPS, see [13]. 
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2.1 MPS formulation of DMRG in finite systems 

Let us start with the MPS formulation of the standard DMRG in finite systems. Throughout 
this paper, we focus on the MPS of open boundary conditions. A spin chain can be bi- 
partitioned into two parts: the system block where the spins are labeled by Sj, and the 
environment block where the spins are labeled by £i. Then we denote an MPS on an L-site 
lattice by the form 

IV^) = Y,A{^ ■ ■ ■ Al^AEl" ■■■E{^\s^--- sps, •••£!) (1) 

where p + q = L and each matrix at each site has dimension m x m with the exception 
that the end matrices A^^ and Ef^ are a row vector and a column vector respectively. 
The A-matrices (E'-matrices) satisfy the left (right)-normalization condition ^Af^Ap = I 

Si 

iYl^i'^i'^ = -^)- It is very useful to express an MPS in a tensor network representation 

(see Schollwock [13]) as in Fig. 1(a). For example, the set of matrices Ap" at the p-th. 
site of the system block corresponds to a tensor (Ap)ap_i,ap with two bond indices ap_i 
and Up and one physical index Sp. The bond index labels the reduced orthonormal basis 
state of the system block at different lengths and the physical index labels the state of the 
spin. Thus the center matrix A represents in fact the wave function of the superblock in 
the reduced orthonormal basis {|apCg)}o,^=i.- - .^i (sec Fig. 1(a)). The essential notion of 
the MPS formulation of the finite-size variant of DMRG is that of local updates; that is, we 
free a matrix, say Eg'' , at a time while keeping the others fixed, update it by A^''^-^ through 
an efficient optimization of the total energy, and shift the center matrix to the right by one 
site. When all matrices in the MPS are updated once, it is called a sweep. Repeat the 
sweep until convergence is reached. 

2.2 E«S»E scheme and the related BMPS representation 

Now consider a special bi-partitioning of the chain as indicated by the notation E»S»E 
where the system block is taken as a consecutive segment of sites around the center of 
the chain and the environment block is taken as two equal separate segments surrounding 
the system block. The superblock is updated at a time by freeing (or adding) two spins 
in between the two sub-systems so that the matrices of the MPS must be associated with 
two distant sites. This special configuration is shown to be particularly adapted to the 
calculation of two-point correlation functions of ID quantum systems or 2D classical lattice 
models [11]. Actually, the MPS for this E«S»E scheme can be formulated as equivalently 
as the MPS for the standard DMRG. By folding the spin chain from the center so that the 
two sub-systems and the two free spins are aligned, every two aligned spins can be regarded 
as a single big spin and the resulting new chain has the same configuration as the more 
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standard S»E scheme. To avoid the notation from getting too involved, it is convenient 
and comprehensive to represent this MPS as a tensor network as in Fig. 1(b) where the 

tensor iAp)'^^__^ ,^^ , now, with two physical indices representing the states of the two spins. 
Alternatively, following McCulloch [5], such an MPS can be simply expressed as 

\i;) = A^---A^AE-...E[ (2) 

In this notation, although the ket basis vectors of the MPS are suppressed and the alphabets 
labeling the state of the aligned spin pairs are simply denoted by two dots, all information 
are preserved. It is easy to understand that, here, the MPS matrices are local state valued 
albeit not explicitly written. 

When a non-Hermitian transfer matrix T is considered, the left and right Perron state 
are generally distinct. Suppose the left and right \ip) Perron state have MPS as follows. 

\^)=B[...B'-prF'-g...F[ 

When we free the two spins associated with the matrices E'^ and F'^, we are actually carrying 
out the local updating by maximizing the quantity {tp\T \ — X {xp\ ip) . Hence, if we impose 
the so-called left- and right-biorthonormal conditions on the MPS matrices 

jy^A, = jy^F^ = I (4) 

then the optimization will be equivalent to an eigenvalue problem of the reduced trans- 
fer matrix T that can be expressed as a tensor network as in Fig. 2 (see Schollwock 
[13]) where the transfer matrix T is expressed as a Matrix-Product-Operator (MPO) and 
ai and gr label the states of the free distant spin pair associated with the matrices Eg 
and Eg. The matrix elements of the reduced transfer matrix can be written explicitly as 

TapaLani^_-„pyyj,C^-i = {(^P<^L(^R^q-i\T \ /3 pC^WnC q -i) ■ ^^c pair of MPS states Eq. (3) 
satisfies the biorthonormal conditions Eq. (4), we refer to the MPS as a dual BMPS. Thus, 
the BTMRG algorithm can be perfectly fitted into the framework of MPS analysis and 
enjoy the same numerical complexity as the Hermitian case in DMRG algorithm. In this 
paper, we will propose two alternative methods for determining the BMPS for the Perron 
state and compare their numerical performance. The details will be described in section 3. 

2.3 Density matrix and canonical form of BMPS representations 

Unlike the Hermitian Hamiltonian in DMRG that there is a unique normalized ground state, 
there are two distinct Perron states in TMRG so that the density operator for the Perron 

state of the system must be taken as p = \(p) (with proper normalization (p) = 1). 
Thus, from Eq. (3), we can readily obtain the reduced density operator for the system and 
environment block as 
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Note that {ap\ Pp) = 6 a/} and Cq) = <^$c view of the biorthonormal conditions. How- 
ever, given a BMPS as in Eq. (3), we can always construct another BMPS by introducing 
an arbitrary invertible (non- unitary) transformation to the matrices E'^ and F'^ (i-e., ap- 
plying a non-unitary basis transformation X to the current biorthonormal bases and 
{\Cq)}) such that Wg = X-^E'g and T'g = X^F^ remain right-biorthonormal. Similar 
results A'p = ApY'^^ and Wp = RpY are obtained by applying a transformation Y to the 
current biorthonormal bases {|ap)} and {|/3p)}. Thus the varied BMPS turns out to be 



I \tlj)=A,---A'-^A'E'-g...E, 

\ \^) = B\...B'-^r'F'-g...F\ ^ 

where the new center matrices become A' = Y^ AX and F' = Y~^FX~^^ and the new 
reduced density operators are simply a similar transform of the old density operators: 'pg = 
Y^^'pgY and = X~^^^X. This implies that if we employ a suitable transformation, we 
can canonize the form of the BMPS representation. 

Given the BMPS in Eq. (3), assume the density matrix pg = FA = YDY~^ is 
diagonalizable, then it can be readily obtained that = F^ A = XDX~^ where X = 
F^Y~^K Since X and Y are unique up to a scaling, we can replace X and Y by XXi 
and YX2 where Ai and A2 are two diagonal matrices. Accordingly, we have new center 
matrices A' = A2-DA1 = /D and T' = A2 ^A^^ = by choosing A2 = .0-^2^-1. Thus 
the canonical form of BMPS reads 

/ \ij) = A,---ApVDE'-g---E\ 



[ \<p) = B,---B-pVDF-g---F, 

and the density matrix bears the canonical form pg = p^ = D. Surprisingly, according to 
our practical simulation, for large enough system size, Pg = YDY~^ always exits and the 
eigenvalues are always real non-negative. 



3 Biorthonormal MPS analysis of iBTMRG in the thermo- 
dynamic limit 

The iDMRG grows the system at a time by adding one or two sites at the center of the 
lattice and produce a fixed point of IMPS that is translationally invariant. By thinking 
each aligned spin pairs as a big spin, our iBTMRG with the E«S«E configuration can be 
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formulated as equivalently as the iDMRG. Throughout this paper, we choose the "single- 
site" iBTMRG algorithm based on the following reasons. First, "two-site" iBTMRG means 
adding four spins at a time where the dimension of the reduced transfer matrix significantly 
increases. Second, the inevitable non-zero truncation of the wave function of the two- 
site iBTMRG often leads to a less- well converged wave function [5]. Third, the single-site 
iBTMRG can achieve zero truncation of the wave function and the problem of being trapped 
in local minimum can be avoided by introducing White's correction [14] to the reduced 
biorthonormal bases. 

In this paper, many formulations are derived from a common basic procedure which we 
will refer to as the Biorthonormalization Procedure. Given two arbitrary bases {|a)}a=i,- ,m 



and {|^)}^=i,...,m , suppose A 



'a=l,-- 



and B = 



are two matrices 



formed with \a) and as their columns respectively. By carrying out the Singular- 
Value-Decomposition (SVD) A^B = USV\ we can readily obtain {AP)\BW) = I where 
P = UE~^/'^ and W = VS~^^'^. This means that, by applying the non-unitary basis 
transformation P and W to the original bases {|a)} and we can obtain a dual 

biorthonormal bases AP and BW. 



3.1 Two alternative iBTMRG algorithms in the thermodynamic limit 

We first propose an iBTMRG algorithm where the main steps are depicted in Fig. 3. The 
iBTMRG starts with an initialization procedure (i.e., n = 1) which selects a minimum length 
of the lattice (depends on the number m of states kept by BTMRG), obtains the respective 
reduced bases of the two sub-systems from the Perron states, and proceeds to biorthonor- 
malize them. The routine in Fig. 3 grows the lattice by adding first two sites to the system 
block and then adding two sites to the environment block, which constitute a period of 
the whole algorithm with respect to a repeated fragment of the final iBMPS. At each time 
when adding new sites to the lattice, we obtain the Perron states (e.g., V ) in the reduced 
biorthornormal bases, carry out an SVD for the Perron states to obtain the reduced bases 
of the enlarged block (e.g., tjj' = 74^^^yl„+i where ^^n+i^n+i = 1)^ ^^d biorthonormalize 

them to build a new biorthonormal bases of the enlarged block (e.g., the right bond states 
of the tensors A^+i and B'^_^_l). If the lattice size grows large enough, we can additionally 
canonize the BMPS. This step is convenient for the purpose of checking the convergence 
of the algorithm. According to the previous description, the biorthonormalization and the 
canonization procedures are simply to apply two successive non-unitary basis transforma- 
tion to the bases obtained from the SVD, we can combine these transformation together and 
denote them by P, W, Q, and R with respect to A, B, E, and F matrices respectively (e.g., 
A'n+i = A]^^^Pn+i and E'^+i = Qn+iE'^+i, and thus An+i = P'^^An+iQ'li). Note that, 
during the SVD in step 2, we keep zero truncation of the wave function when obtaining the 
reduced basis of the enlarged block. The resulting BMPS read 
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\ip) = B\---B-^rnF-^---F, 



(8) 



The above algorithm builds the biorthonormalized bases from the reduced bases derived 
from the SVD that are based on another biorthonormal bases. One may notice that there 
exist actually a lot of possibilities to build the biorthonormal bases. Nevertheless, the 
process of SVD is critical. We can also jump the SVD and directly biorthonormalize t/'' 
and ip'. But such a numerical procedure will become unstable very quickly. Here, we 
provide another algorithm which is very different from the above algorithm in nature. The 
alternative algorithm is depicted in Fig. 3 with the steps 2-3 enclosed in a rectangle being 
replaced by the routine depicted in Fig. 4. To the contrary, this algorithm constructs the 
biorthonormalized bases from the reduced bases derived from the SVD that are based on 
two orthonormal bases. The first major change lies in step 2 where before performing SVD 
the Perron states must be restored to be in terms of their own orthonormal bases, e.g., -(/;'• = 
Pnip'Qn- The second change lies in step 3 where before performing the biorthonormalization 
procedure the reduced bases derived from the SVD must be transformed back to be in 
terms of the previous biorthonormal bases, e.g., = P~^A]^_^_i. As it turns out, two 

remarkable properties emerge, e.g., = P~-'^A^_|_iPn+i and A'^ ■ ■ ■ A,^ = A[ ■ ■ ■ A^Pn^i. 

Similar relation holds for the B, E, and F matrices. Accordingly, the BMPS can also be 
expressed in terms of two series of reduced orthonormal bases. 



where An = PnAnQn and = WnFnRn- For convenience, in this paper, the former 
algorithm will be referred to as iBTMRG A and the latter iBTMRG B. 

3.2 Efficient wave function prediction and the fixed point of the iBMPS 

In both algorithms, the iBTMRG uses a local eigensolver to obtain the Perron vectors. 
A good initial guess for the eigenvector will significantly improve the performance of the 
eigensolver. In [5] , McCuUoch has developed a wave function transformation in iDMRG for 

the prediction of the wave function as the lattice size is increased. Such a transformation 
can be effectively translated to our iBTMRG. Therefore, at the end of the first half-period 
of the process in Fig. 3, the prediction of the next Perron vectors will be 



where and ■i/'' are the current canonical center matrix and the current Perron vector. 

Similar results go to the Perron vectors prediction at the end of the second half-period in 




(9) 




(10) 
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Fig. 3. 



Atrial = A^^^E-^^,D-'/^D]II^ = ^■■D^^I''d]11 



I • p Z7- n-V2nl/2 ■n-V2nl/2^ 

y Atrial - n+lJ'n+l^n - 'P 

1/2 

Note that, for finite fixed n, the center matrix Dn may not be identical for the two half- 

1/2 

periods in Fig. 3 so that the center matrices in Eqs. (10-11) should use their own Dn ■ 

1/2 — — 

When the size of the lattice is small, the Dn can be replace by An and Fn- 

At the fixed point of the iBTMRG, the translational invariant iBMPS should bear the 
canonical form 

i m = --- (z?;!f A;Dy'^:,)(D;yi'A:,z?y%) . . . 
\ = • • • (z?;yiX^yX)p;-(X^n/X) • • • 

where we take n to be the iteration step when the convergence criterion is met. At the 
fixed point, the density matrix would be invariant no matter where the bond is located. 
This leads to the fixed point criterion Yl^n^nB'n = D^-i- Since Tr{Dn) = 1 (with proper 

normalization (^'| ip') = 1 ), one way of measuring the closeness of the two density matrices 
Dn-i and D!^ = J2A'^DnB'^ is given by the Kullback-Leibler divergence 

I)KL = l:^n-l(^)log(%;^) (13) 

In this paper, Eq. (13) is utilized to check the convergence of our iBTMRG. 

3.3 Normalization of the iBMPS 

The fixed point criterion Y^A^DnB'^ = D^-i is equivalent to A^d]J'^){D~]/^ B'^Dn"^)^ = 

— 1/2 — ■ 1/2 —1/2 — ' 1/2 

/ which means the two matrices ^n^n and D^_^ B'^Dn satisfy the right-biorthonormal 

condition. Similarly, Dn"^ E'^D^]/^ and D^"^ F'^D^]/^ must satisfy the left-biorthonormal 
condition. Thus, the overlap of the iBMPS in Eq. (12) must be exactly (V'| = 1 in the 
fixed point. The normalization = 1 is important when we want to apply the iBMPS 

to the calculation of expectation value or correlation functions in the thermodynamic limit. 
However, in practical simulations, the finite number m of states kept by BTMRG and finite 
iterations can only reach an approximation of the fixed point so that the iBMPS will not 
be exactly normalized. In [15], Orus and Vidal have developed a method for the orthonor- 
malization of an IMPS in iTEBD, which can also be effectively translated to our iBTMRG. 
Define two transfer operator Tr and Tl as shown in Fig. 5. In order to achieve the right- 
biorthonormal condition, assume Qr is the right dominant eigenvector of Tr and an SVD 
allows to decompose Qr = UrVr. We now perform two similarity transformations to the 
unit cells in Eq. (12) which lead to 
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Similar arguments apply to and its left dominant eigenvector fi^ = UlVl and lead to 
the left-biorthonormal iBMPS 



(15) 



[ W) = --- (vis;Dy2F;z);l(V^-i)(yii?;i?y'F;i?;:yiV^-i 

Accordingly, update the matrices: A^ -ir- U^A^ , ^lB'h , E'^ ^ ^n^R ■> ^ ^n^R 

1/2 1/2 

, and ^ VlD^_^Ur. Then the lattice can be separated into two parts at any bond 

location where the iBMPS can be expressed as similar to as Eq. (7) 

f |v> = • • • (A;,z)y2^;z);lf )i)y_\(z);:!f A^z^y^^j . . . 

\ 1^) = . . . (s^^y^F^i^^if )z)y_\(D;if s^^y^Fj ... 

4 Example: 2D classical systems 

In this section, we test our iBTMRG A and iBTMRG B algorithms on a 2D classical 

Ising model. The target matrix is IlJ^ where 77y is the non-Hermitian fundamental transfer 
matrix of the general local energy function (LEF)-parameterized Markov random field on 
an infinitely-long vertical twisted cylindrical lattice with peripheral length N [10, 16] (The 
Ising model is just a special case of a Markov random field). This matrix is intimately 
related to a 2D Markov additive process and enjoys a very special SVD structure and many 
fascinating properties [16]. 

4.1 iBTMRG in finite systems 

Our iBTMRG algorithms can also be applied to the partition function calculation of classical 
models in finite lattices. The method to build the reduced biorthonormal bases in iBTMRG 
A is essentially the same as the previous BTMRG in [11] which we will refer to as iBTMRG 
C for convenience. However, the system-growing strategy of iBTMRG C is different from 
that of the single-site iBTMRG. Upon the same E»S»E configuration, the previous BTMRG 
grows the system by adding two sites at each iteration to the system block while keeping the 
environment block fixed and small. To our knowledge, both proposed single-site iBTMRG 
algorithms had never been tested before. Here, we will compare the performance of our 
iBTMRG algorithms and the previous one by computing the free energy of the Ising model 
on a finite 2D lattice. Figure 6 shows the error of the free energy (i.e., the logarithm of the 
Perron root of the transfer matrix Uq) of an anisotropic Ising model [Jx = —Jy where Jx 
and Jy represent the interaction between horizontally and vertically neighboring spins) at 
the critical temperature = 2/log(l -|- \/2) for system size A'' = 160 and various number 
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of states m kept by iBTMRG. We also compare the performance of the finite-size variant of 
BTMRG (fBTMRG) for the three algorithms where the Perron states from the iBTMRG 
were used as the initial BMPS states and were further optimized variationally. Prom Fig. 
6, we can see that iBTMRG A and B exhibit almost identical performance which were 
further slightly improved by their variational finite-size variants. For fBTMRGs, algorithms 
A and B give rise to very close results and only outperform algorithm C a bit. Although 
iBTMRG C has the merit of high efficiency, it is interesting to note that its accuracy is 
far worse and remain constant for all values of m. This is because, in the system-growing 
stage, the size of the environment block is kept fixed and small so that the entanglements 
between the two subsystems remain fixed and small irrespective of the reduced basis size 
m. 

4.2 iBTMRG in the thermodynamic hmit 

In iBTMRG, finding the Perron states by an iterative eigensolver is the most time-consuming 
part so that the overlap between the initial wave function and the variational optimum 
i'^'triall '^'opt) ^ill dominate the performance of an iBTMRG algorithm. Figure 7 shows the 
fidelity 1 — {''P'triall '^'opt) (s^® McCulloch [5]) between the Perron states predictors and the 
optimal ones for the above anisotropic Ising model at criticality with m = 40 states kept 
in the reduced basis. Before the peripheral size of the system grows larger than N = 50, 
the center matrices in Eqs. (10-11) must be replaced by the un-canonized center matrix 
An and since the density matrices cannot be diagonalized to real positive eigenvalues 
until N = 50. From Fig. 7, the un-canonized center matrices seem inappropriate for the 
usage of wave function transformation. Once the canonical center matrices are utilized, 
the fidelity quickly drops down to 10~® and continues to decay in a nearly constant rate. 
Moreover, both predictions for algorithms A and B have very close effectiveness and the 
right Perron state transformation appears to prevail with respect to the left Perron state 
transformation. Another important issue concerns the convergence of the translationally 
invariant fixed point of the iBMPS. In this paper, we use the Kullback-Leibler divergence 
to monitor the convergence of the algorithm. The inset of Fig. 7 shows the convergence 
to the fixed point of the iBMPS with respect to the number of iterations. In view of the 
degeneracy of the Perron root and the long-range spin correlation at the criticality, the 
convergence appears to be quite slow. However, in off-critical regions, the convergence can 
be as fast as reaching 10~^^ as the size exceeding N = 200. 

Once we have the translationally invariant iBMPS (after normalization {tp] (p) = 1), 
the calculation of the expectation value or the two-point correlator in the thermodynamic 
limit can be very efficient. For illustration, we take an isotropic (J^ = Jy = J) Ising 
model as an example and calculate its magnetization as a function of the temperature 
and its spin-spin correlation function at the critical temperature. Now consider Eq. ^16), 
if we free the unit cells V'^;^/^'"^"^ = ^^.^(D^lf :4^^"«DyX^"«) and <^;;^/^'"«"« = 



11 



D]l^'^{D^]/^B"^^'^'''Dn'^F'^^^'^^) then we have Perron states in the thermodynamic hmit 
expressed as ^^^^^-''^^^^^ and (^^^^^^'^^^^ with respect to a dual biorthonormal bases. Note 
that the spins (Tl(T^ and CFRa'j^ are two neighboring spins pairs located on the left and right 
side of the system block, and the spins cflc^r and cr^cr^ are two (infinitely) far distant 
aligned spins. Therefore, the magnetization in the thermodynamic limit can be calculated 

by 

{-r)=y.-^ e e (17) 

For an infinite 2D isotropic Ising model, the exact solution for the magnetization is ^1 — (sinh(2 J/ZcbT)) 
[17]. Here, we compute the magnetization by Eq. (17) and the exact formula as in Fig. 8 
where the inset shows the relative error between the numerical result and the exact one. 
Similarly, the spin-spin correlation between two infinitely far distant spins in the thermo- 
dynamic limit can be calculated by 

{oLan) = E ^^^^ E E V'::/^"^^'^^;^"'^^'^'^ (18) 

At the critical temperature, the exact spin-spin correlation function scales as G{r) oc 
j,-i/4 ^-y _|_ (9(7--2-j^ ^-|^gj_ Ideally, Eq. (18) will be zero but in practice it can only be seen 
as the correlation of two spins separated by a distance r equal to half of the system size 
N where the fixed point criterion is met. This implies that we can regard Eq. (18) as the 
correlation function in the thermodynamic limit as r = N/2. Figure 9 shows the scaling 
behavior of the spin-spin correlator compared with the exact solution. It is worthy of not- 
ing that, in the low-temperature region nearby the criticality, our iBTMRG behaves quite 
prone to getting stuck. Fortunately, by introducing White's correction [14] to the reduced 
biorthonormalizcd bases, the algorithm exhibits much improved efficiency and convergence 
with the correction weight around lO^'^-lO^^. In [15], Oriis and Vidal have conducted the 
same calculation as in Figs. 8-9 by using the iTEBD algorithm for the same Ising model 
(see Figs. 11-12 in [15]). For the magnetization, although the iTEBD has achieved a better 
accuracy than iBTMRG, they are not on the same conditions. In addition to the main 
difference already mentioned in the introduction, another difference is that the transfer ma- 
trix they have dealt with is Hcrmitian which is much numerically well-conditioned than the 
non-Hermitian case. More importantly, for the two-point correlation calculation, instead of 
evaluating a very long tensor network as in [15], we use a very efficient formula Eq. (18) 
and obtain a result that is comparable with respect to the iTEBD. 
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5 Conclusions 



In this paper, we give a thorough Biorthonormal Matrix-Product-State (BMPS) analysis 
of the Transfer-Matrix Renormahzation-Group (TMRG) for non-Hermitian matrices and 
propose a "single-site" infinite-size Biorthonormal TMRG (iBTMRG) algorithm to directly 
obtain a dual infinite-BMPS (iBMPS) which are translationally invariant in the thermody- 
namic limit. Different from the standard DMRG that the MPS was established on a series 
of reduced orthonormal bases, the BMPS are built on a dual series of reduced biorthonor- 
mal bases for the left and right Perron states of a non-Hermitian matrix. We propose two 
alternative methods for the construction of the dual biorthonormal bases and compare their 
numerical performance in both finite and infinite systems. We also develop an efficient wave 
function transformation of the iBTMRG, an analogy of McCulloch [5] in the infinite-DMRG, 
to predict the wave function as the lattice size is increased. The resulting translationally 
invariant iBMPS not only allows for investigating bulk properties of a strongly correlated 
system in the thermodynamic limit without the boundary effect but also allows for the 
evaluation of expectation and two-point correlation function of the system very efficiently. 
For illustration, we calculate the magnetization as a function of the temperature and the 
critical spin-spin correlation function in the thermodynamic limit for a 2D classical Ising 
model. 
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FIG. 1. (a) The tensor network representation of the MPS in the standard S«E configuration, (b) The tensor 
network representation of the MPS in our special E • S • E configuration. 




FIG. 2. The tensor network representation of the reduced transfer matrix T associated with the dual BMPS of 
BTMRG upon the E«S«E configuration. 
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Step 3 . I Perform Biorthonormalization 
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FIG. 3. The main steps of iBTMRG A, which grows the lattice by adding first two sites to the system block and 
then adding two sites to the environment block, constitute a period of the algorithm with respect to a 
repeated fragment of the final iBMPS. 
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FIG. 4. The main steps of iBTMRG B is the same as iBTMRG A only with the steps 2-3 enclosed in a rectangle 
in Fig. 3 being replaced by the routine depicted in this figure. 
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FIG. 5. To achieve the right (left)-biorthonormalization condition, we must perform two similarity 
transformations to the unit cells of the iBMPS where = ^^r {^l the right (left) dominant 

eigenmatrix of the transfer operator T^, ( ). 




FIG. 6. The comparison of the infinite-size (iBTMRG) and finite-size (fBTMRG) variants of three BTMRG 
algorithms by plotting the error of the free energy of an anisotropic Ising model with lattice size A'=160 at 
the critical temperature T^, =2/log(l + V2). 




lattice size 

FIG. 7. The fidelity l-^v|/|'j.j^[|v|/^p,^ between the Perron states predictors and the optimal ones for the anisotropic 

Ising model at criticality with m=40 states kept in the reduced basis. The inset shows the convergence of the 
Kullback-Leibler divergence to the fixed point of the iBMPS with respect to the number of iterations. 
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FIG. 8. The magnetization as a function of the temperature of an isotropic Ising model in the thermodynamic 
limit with m=40. The inset shows the relative error of the magnetization with respect to the exact solution. 
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FIG. 9. The scaling behavior of the spin-spin correlator of an isotropic Ising model at the critical temperature 
obtained from Eq. (18) which is compared with the exact solution. 



